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Abstract. We investigate models of dark energy with purely kinetic multiple k- 
essence sources that allow for the crossing of the phantom divide line, without violating 
the conditions of stability. It is known that with more than one kinetic k-field one 
can possibly construct dark energy models whose equation of state parameter w x 
crosses —1 (the phantom barrier) at recent red-shifts, as indicated by the Supernova la 
and other observational probes. However, such models may suffer from cosmological 
instabilities, as the effective speed of propagation c x of the dark energy density 
perturbations may become imaginary while the w x = — 1 barrier is crossed. Working 
out the expression for c x we show that multiple kinetic k-essence fields do indeed lead 
to a w x = — 1 crossing dark energy model, satisfying the stability criterion c? > as 
well as the condition c x < 1 (in natural units), which implies that the dark energy 
is not super-luminal. As a specific example, we construct a phantom barrier crossing 
model involving three k- fields for which c x is a constant, lying between and 1. The 
model fits well with the latest Supernova la Union data, and the best fit shows that 
w x crosses —1 at red-shift z ~ 0.2, whereas the dark energy density nearly tracks the 
matter density at higher red-shifts. 
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1. Introduction 

The late-time accelerated expansion of the universe is now well established through 
the observations of the luminosity distance - red-shift relation for type la Supernovae 
(SNIa) [1, 2, 3, 4, 5, 6, 7, 8]. Measurements of the cosmic microwave background (CMB) 
temperature fluctuations by the Wilkinson Microwave Anisotropy Probe (WMAP) 
[9, 10] as well as the large scale red-shift data from the Sloan Digital Sky Survey 
(SDSS) [11] further indicate that the universe is very nearly spatially flat and consists 
of about 70% of dark energy (DE) which drives the cosmic acceleration. The value of 
the DE equation of state (EoS) parameter, i.e., the ratio of pressure to energy density 
of DE, w x = p x /p x , should be less than —1/3 to maintain this acceleration. Although 
the simplest candidate for DE is a cosmological constant (A), for which w x = — 1, 
there are serious theoretical problems associated with it [12]. A needs to be extremely 
fine tuned so that the DE at present is small compared to the Planck scale. A fine 
tuning is also required to keep the (A-) dark energy density p x comparable with the 
present critical density p 0c , while the acceleration begins only in the recent past - 
the coincidence problem. From the observational point of view, although there are 
indications that w x is very close to —1 at present, combined studies of SNIa, CMB 
and baryon acoustic oscillation (BAO) peaks do actually point toward a dynamical DE 
with EoS parameter w x slowly varying with red-shift z. More specifically, the so-called 
ACDM model, which consists of A and cold dark matter as dominant components, 
although consistent with the low- ,2 observational data, does not produce sufficiently 
good fits with the data at relatively higher z (> 1). Accordingly, there have been a large 
number of model-independent and model-dependent DE parameterizations explored in 
the literature [13, 14, 15, 16, 17]. Extensive studies reveal that in parameterizing the 
DE to produce good fits with the observational data, it is always desirable to have no 
inherent restriction imposed on w x from the underlying theory. For instance, it is not 
desirable to have w x always > — 1, which is required for the theoretical consistency of 
many popular DE models, such as quintessence, tachyon, dilaton, etc. For the reviews 
on various dynamical aspects of DE, see [18]. 

Within the scope of minimally coupled canonical scalar field DE models, such as 
quintessence [19], the crossing from w x > —1 to w x < —1 can never happen and there 
is a need of more than one scalar field, at least one of which is of phantom nature, i.e., 
carries a wrong sign in front of the kinetic term [20]. For this reason the w x = — 1 
line is commonly known as the phantom divide line (PDL). Phantom fields violate the 
dominant energy condition, and therefore could give rise to classical instabilities [21]. 
Moreover, since the phantom energy density is unbounded from below, the quantum 
mechanical stability is in jeopardy. The vacuum is unstable against the decay into ghosts 
and positive energy particles, which couple to gravitons producing further decay into 
cosmic gamma rays [22]. In an expanding Friedmann universe, the extremely repulsive 
(anti-gravitating) nature of the phantom fields pose more problems as the phantom 
energy density generically tends to infinity in finite time in future leading to the so- 
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called "Big Rip" singularity [23] , as well as some other mild singularities. For a general 
categorization of these singularities see for instance [24]. 

Apart from the coupled quintessence and phantom models, commonly known 
as the quintom models [25], there are various other models in which the phantom 
barrier crossing is shown to be achieved in a consistent way. Notable among these are 
the scalar-tensor models [26], brane- world models [27], string-inspired dilatonic ghost 
condensate models [28], quantum-corrected Klein-Gordon models with quartic potential 
[29], coupled DE models [30], modified gravity models [31], H-essence (complex scalar) 
models [32], etc. In this work, we focus on the so-called k-essence models [33, 34], which 
incorporate a natural generalization of quintessence via a Lagrangian with non-linear 
dependence on the kinetic term(s) for the scalar field(s). For such models w x is not 
constrained to be > —1 and may in principle cross the PDL. However, it has been argued 
by Vikman [35] that in the case of a single field k-essence model, with a Lagrangian of 
the form P(<f), X), where X = — |[V0] 2 is the kinetic term of the scalar field 0, the PDL 
crossing either leads to instabilities against cosmological perturbations or is realized by 
a discrete set of phase space trajectories^. For a multiple k-essence model, involving 
more than one k-fields 7 (/ = 0, . . . , TV — 1), with kinetic terms X 1 = — |[V0 7 ] 2 , the 
PDL crossing is shown to be possible even when the Lagrangian is purely a function of 
the kinetic terms, P = ^PjiXj) [38]. However, the stability of such a model against 
cosmological perturbations have not been explored in the literature. 

A major difficulty in dealing with the equations of scalar-type cosmological 
perturbations for multiple scalar fields is that these equations are extremely coupled 
in the curvature and isocurvature modes, even when the scalar fields are only minimally 
coupled to gravity [39, 40, 41, 42]. For multiple k-essence, the perturbation equations 
are complicated further, however there are some specific cases in which the expression 
for the speed of propagation c x of the curvature perturbations could be obtained 
and the stability analysis could be performed [43, 44]. For example, when the k- 
essence Lagrangian is P($,X), where <l> = {4> ,4> 1 , ■ ■ •) and X = — i^ /J ($)V M / V /1 J 
[43], the propagation speed is simply given by the well-known result of Garriga and 
Mukhanov [45], i.e., c 2 x = ^/-§x~- On the other hand, in the multi-field Dirac- 
Born-Infeld (DBI) inflationary scenario, the propagation speed for the Lagrangian 
P = — / _1 ($)a/1 — 2X~JJW) — V(<f>), is given in a homogeneous background by 
c x = y/1 - 2XJW) [44]. However, none of these expressions for c x will hold for the 
most general Lagrangian P = ^P I {4> I -,X I ), corresponding to the N non- interacting 
fields $ = (0 O , (f) 1 , ■ • •). A comprehensive study of the curvature perturbation equations 
reveals that one can, in principle, obtain the expression of c x for the multi k-field 
Lagrangian P = ^2 P } ((j> } , X } ) [46], leaving aside the isocurvature modes which couple 
weakly with the curvature mode in the large-scale limit [40, 41]. 

In this paper we first investigate whether the multiple kinetic k-essence DE models 

| There are however some counter-arguments and counter-examples set in [36, 37] and it is shown that 
for a generalized single field k-essence Lagrangian of the form P(cf>,X) — f(<fi)X — V(4>) stable PDL 
crossing is indeed possible for some specific functional form of f(<j)) and V(cf)). 
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that allow the PDL crossing is cosmologically stable, i.e., the speed of propagation c x of 
the DE density perturbations satisfies the condition c 2 x > 0. It is also relevant to check 
whether c x < 1 (in natural units), so that the DE perturbations do not travel faster 
than light [47]. However, the requirement c x < 1 is not potentially stringent, because 
it has been argued that causal paradoxes do not necessarily arise even when c x > 1 
[48]. We find that in most circumstances (including the model described in [38]) where 
c x dynamically varies with time, the PDL crossing is associated with a flip of sign of 
c 2 x , thus giving rise to violent instability. An appropriate choice for c 2 x in a stable DE 
model is therefore a constant whose value lies between and 1. Imposing this condition 
on c? we work out the constraint on the field configuration. As a specific example we 
try to construct a stable PDL crossing DE model using this constraint and with the 
least number of model parameters possible, so as to determine these parameters with 
less uncertainties by fitting with the observational data. While within the scope of 
two k-fields the stability constraint leads to complicated field solutions, we construct 
a reasonably simpler model involving three k-fields for which stable PDL crossing is 
possible. 

We constrain the parameters of the model using the most recent data-set released 
by the Supernova Cosmology Project (SCP) team [7], viz., the "Union" compilation, 
which consists of 414 SNIa reduced to 307 usable data points after various selection 
cuts. Marginalizing over the Hubble constant H , we find that the best fit value of the 
present matter density parameter Oo m ~ 0.265, whereas the best fit DE EoS parameter 
w x crosses the PDL at red-shift z = z c ~ 0.2, fairly in agreement with z c found in many 
model-independent DE parameterizations [49]. The minimized value of the x 2 statistic 
is found to be about 310.4, which shows an improvement of A% 2 ~ —2.7 over the 
minimized \ 2 f°r the ACDM model found with the Union data-set [8]. The minimized 
X 2 per degree of freedom for the model discussed in this paper is 1.011, which is also 
marginally better than that found in [6] for the ACDM model using the older SNIa 
data-set complied by Riess et al [4]. Additionally, we see that the best fit value of the 
DE density almost tracks the matter density over sufficiently longer red-shift range in 
the past until becoming dominant after z ~ 0.4. Thus the model provides a possible 
resolution of the coincidence problem. Extrapolation of the best fit w x to future epochs 
shows a tendency of w x to turn back towards —1 after reaching a minimum value. Thus 
the phantom regime appears to be transient and future singularities, such as the Big 
Rip, may be avoided. 

The present model also differs from the usual quintom model [25] implemented 
using two (or possibly more) scalar fields, (at least) one of which has a negative definite 
kinetic energy density (i.e., phantom), and as such the EoS parameter is less than — 1, 
whereas the other field(s) has(have) positive definite kinetic energy (quintessence-type) 
and EoS parameter (s) greater than —1. The negative energy field comes to dominate 
at late times and is responsible for the PDL crossing. In the three-field k-essence model 
which we study, two of the k-fields have positive definite energy density, whereas the 
energy density of the third k-field is negative for a sufficiently longer time in the past 
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and present, but may become positive in distant future, as found by extrapolating the 
best fit model parameters to future (negative) red-shifts. In that sense the third k-field 
may not be regarded as equivalent to an eternal phantom, and quantum mechanical 
instabilities may be avoided. None of the fields have EoS parameter less than — 1, 
unlike quintom models. Moreover, the negative energy k-field is not responsible for 
the super-acceleration (w x < —1), it is a positive energy k-field that causes the PDL 
crossing, again in contrast with what happens in quintom models. 

This paper is organized as follows: in sec. 2 we describe the general framework of 
multiple k-essence DE models assuming that there are no mutual interactions between 
the individual k-fields. In sec. 3 we consider the simplistic case of purely kinetic multiple 
k-essence and show that only for more than one dynamical k-fields the PDL crossing 
is possible, although the field configuration is severely constrained by the stability 
criterion. We work out this constraint and choose the value of the speed of propagation 
c x of the DE density perturbations to be a constant, in order to have a stable PDL 
crossing multiple kinetic k-essence DE model. In sec. 4 we construct a specific multiple 
k-field model for which c x is constant, and w x crosses the PDL in recent past. In sec. 5 
we fit this model with the 307 Union SNIa data [7], to reconstruct phenomenologically 
the multiple k-essence Lagrangian. We conclude with a summary and open questions 
in sec. 6. In the Appendix, we work out the equations of scalar-type cosmological 
perturbations for a multiple k-essence field configuration so as to find the expression for 
the propagation speed c x . 

2. General Formalism 

Let us consider the following action in (3 + 1) dimensions: 

N-l 

P(*,X) = ^P / (0 / ,X I ), (1) 

1=0 

where k 2 = %nG is the gravitational coupling constant, C m is the Lagrangian density for 
perfect fluid matter, <3> = (0 O , • • • , (p N _ 1 ) are a set of N mutually non-interacting scalar 
fields, and X = (X , • • • , X N _ 1 ) where X 7 stands for the kinetic term of the scalar field 

X, = -Ig^dnfad^ , (no sum over / = 0, • • • , iV-1) . (2) 

Pj^jjXj) is the k-essence Lagrangian density for the field 

We assume the space-time structure to be described by a spatially flat Robertson- 
Walker (RW) line-element 

ds 2 = -dt 2 + a 2 (t) [dr 2 + r 2 (d$ 2 + sin 2 $dip 2 )] (3) 

where a(t) is the scale factor which is normalized to unity at the present epoch t — t . 
In such a space-time, the kinetic term reduces to X l = <j) 2 /2, where the over-dot stands 



S = d A x 
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for d/dt. The field equations are the Friedmann equations given by 

H 2 = - = — (prn + P x ) , H= = -— [ Pm + p x +p rn + P x ] , (4) 



3\r iti> ' r jt / 7 o o 

a a 2 2 

where H = a/a is the Hubble expansion rate, p m and p m are the perfect fluid energy 
density and pressure respectively, p x and p x are the k-essence energy density and 
pressure respectively, given by 

N-1 N-1 /op \ N-1 

Px = = E( 2X ^ - P A > Px = ( 5 ) 

Considering the perfect fluid matter in the form of pressure-less dust {p m = 0), the field 
equations integrate to give the matter energy density as p m = po m /a 3 , where p 0m is the 
present value of p m . From Eq.(4) one obtains the expression for the normalized Hubble 
rate: 



n l V 2 

Px _|_ "0m 

Poc a 3 . 



(6) 



where if is the present value of the Hubble parameter H and p 0c = 3Hq/k 2 , Q 0m = 
Pom/ Poc are the present values of the critical density and matter density parameter 
respectively. 

The equations of motion of the scalar fields <& = (0 O , • • • , Ar _ 1 ) are given as 

a W, (forall/) - (7) 

and the mutually non-interacting individual field components satisfy the continuity 
equation: 

Pl = 3H(1 +w I )P I (for all /) , (8) 

where Wj = P I jp I are the EoS parameters for the individual fields 7 . The EoS 
parameter w x for the k-essence DE is thus given by 

N-1 

Px_ = WiPi 

Using Eqs. (4) and (6) one can express the total EOS parameter w in terms of w x : 

w = n = w x 1 - • (io) 

Pm + Px \ H 2 a 3 J 

The deceleration parameter q is given by 

a l + 3w 

q= = —2~- 

For an accelerated expansion of the universe (q < 0) the total EOS parameter w 
should be less than —1/3, and as can be observed from Eq. (10), w x is further less than 
— 1/3. Although the observations favour the value of w x to be very close to that for a 
cosmological constant (i.e., —1) at present and in recent past, extensive studies show 
that better fits with observational data require a dynamical DE with w x substantially 



w x = — = > — = "I" • (9) 

X Px £j Px 3Hp x 
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greater than —1 at relatively higher red-shifts (z = a -1 — 1 > 1). Moreover, for consistent 
statistical data fits it is not desirable to always have a theoretical restriction such as 
w x > — 1. In other words, the dynamical DE models that allow the PDL crossing are 
most suited for analyzing the observational data so as to determine the cosmic evolution 
history over the entire red-shift range that is probed. While such a crossing in minimally 
coupled scalar field DE models is associated with violent cosmological instabilities, one 
has to ensure that the instabilities never occur in the multiple k-essence DE models. 
For such models the expression of the square of the speed of propagation of curvature 
perturbations is found to be (see the Appendix for detailed derivation) 



c l = (Px + P> 



'N-l 

E 



p_ 1 ±p 1 

c 2 



-1 



~ 2 dPJdX, 

where c < = i^k;' (12) 



./=o 

In the case of a single k-field (J = 0, say) the above Eq. (12) reduces to the well-known 
result c 2 x = c 2 = / §y~ [45]. The value of c 2 x must be non-negative, for the stability of 
a model, and always be less than unity (in natural units) so that the DE perturbations 
do not travel faster than light. 

In the next section we examine the plausibility of a stable PDL crossing in the 
multiple k-essence scenario, considering for simplicity the k-field Lagrangian to be purely 
kinetic in nature. 



3. Purely kinetic multiple k-essence 

When the k-essence Lagrangian is purely a function of the kinetic terms X = 
(X , • ■ • , X N _ 1 ), i.e., Pj = P I (X I ), (I = 0---,iV— 1), the equations of motion (7) 
for the scalar fields $ = (0 O , • • • , Ar _ 1 ) can be integrated at once to give 
dP, 

constant (kj) , for all /. (13) 



a 



3 i u± I 



1 dX 



i 

'2 



However, remembering that X l = <pj2 in Robertson- Walker space-time, one may 
observe that the expressions (5), (9) and (12), for p x , w x and c 2 respectively, remain 
invariant under a field redefinition 7 — ■> k 1 $ l . Hence the constants k t can effectively be 
set to unity, for all J, and the above equation can be rewritten as 

, Ap 1 

Now, using the definitions (5) of p x and p J , we can express the DE EoS parameter 
w x , Eq. (9), and the square of the propagation speed of perturbations c 2 x , Eq. (12), as 

w = _ 1 + — Y"2X c 2 = J2f=o 2X i dP i/ dX i (15) 

p x j^ 1 dXj - Y.to[^ I dP I /dX I +AX^P I /dX^]' 

For a transition from w x > —1 to w x < —1, the quantity s ^2X [ dP I / dX t needs to 
vanish at the PDL cross-over scale factor a = a c . 

In the case of a single scalar field <pj = 4> , say, and P r (X J ) = P Q (X ), the PDL 
crossing implies that X (dP /dX )\a =ac = 0, i.e., either X (a = a c ) or dP /dX \a =ac or 
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both must vanish. This is clearly not possible for any finite value of a since from Eq.(14) 
we have ^/W {dP /dX ) = a~ 3 . 

For more than one k-fields the PDL crossing is in principle possible whenever the 
quantity ^2X I {dP I /dX I ) = at a = a c , i.e., Y A/2X 7 | a=(lc = 0, by Eq. (14). However, 
from Eq. (15) we find that the change of sign of Yl ZX^dPjdX^ as the w x = —1 line 
is crossed, implies a change of sign of sign of c 2 x as well, unless the numerator and 
denominator of c 2 x exactly cancel each other to give a value which never flips sign. For 
example, if P z is a function of some numerical power of X n i.e., P r = X" 1 (for all J), as 
has been considered in ref. [38], there is always a flip of sign of c 2 x associated with the 
PDL crossing. Therefore unless one chooses c 2 x to be such a function of a which does 
not change sign for any value of a, instabilities are bound to be there even in multiple k- 
essence DE models. Moreover, we must have c 2 < 1 in order that the DE perturbations 
do not travel faster than light. The simplest possible choice of c 2 x , appropriate for a 
stable and non-super-luminal DE model is 

c 2 = n (constant) , < n < 1 . (16) 

Of course, the limiting case n = implies a non-propagating DE (for example, due to a 
cosmological constant). For dynamic DE models we have < n < 1. 
Let us now denote 

y, = \JlX l = 4> T in RW space-time . (17) 

Eq. (14) reduces to 



-r 1 - = -7 for all I , 

ay j a 6 

and the DE energy density takes the form 

JV-l 

Px = P i ' Pi 



1=0 



a 3 1 



From these equations we find 



a p'M) 



yA a ) 



for all I , 



(18) 



(19) 



(20) 



where prime {'} denotes d/da. Eqs. (15) for w x and c 2 reduce to 



-1 



a 3 p x 



ay{a) 



N-l 



1=0 



N-l 



where y(a) = ^y^a) . (21) 



We rewrite the expression for c 2 in the form 



y 

3 c x77 + a ^( fl ) = 



77,(0), 



(22) 



where the arbitrary functions 77,(0) satisfy the relation Y r li( a ) = 0- With the choice 
c 2 = n (constant), we therefore have 



Vi 



1 + 



ayj(a) 



3n 
a 



N-l 



^77,(0) =0. 



(23) 



1=0 
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While from Eq. (21) we see that the PDL crossing requires y = YlVi — at 
stable k- field configuration is the one that solves the above equation (23). 
In the next section we construct a cosmologically stable PDL crossing DE model by 
suitably choosing a specific field configuration and a reasonably simplified ansatz to 
solve for y, w x and p x . We fit the solutions to the latest SNIa Union data [7] in the 
subsequent section, show the plausible tracking of the dark energy density to the matter 
density in the past, avoidance of future singularities, and finally reconstruct the multiple 
kinetic k-essence Lagrangian. 



4. Cosmologically stable PDL crossing multiple kinetic k-essence model 

4-1. Two k- field configuration 

Let us first consider a field configuration involving only two scalar fields: 

7 = 0,1 ; $ = (0 ()! i ) ; rj (a) = -7^ (a) = 77(a), say. (24) 
From Eq. (23) we have 

<(l + ^-)= *(l-l-) = -*i. (25) 

Let us assume a power-law solution for y 1 as 

y 1 = a 1 a s , a 1 , s — constants . (26) 

Plugging this in Eq. (20) and integrating, and using Eq. (19), we get the expressions 
for the energy density and the pressure corresponding to the field <f> 1 as 

« 1 ° S " 3 + ^ ( s ^ 3 ) P = { ~ fe) « 1 qS " 3 " ^ 3 ) (27) 

11 y -3a 1 lna + A 1 (s = 3) ' 1 | a x a s ~ z - 3a, In a - A, (s = 3) 1 ; 

where X 1 is a constant. 

Substituting Eq. (26) in Eq. (25) and eliminating 77(a) we obtain 

V 3nJ y l + yjy, 1 ' 

Integrating this we get the algebraic equation 

y + Vl (1 - Cyl +S ^ n ) = , (29) 

where C is an integration constant. 

We look for the simplest possible solutions of the above equation (29), so that the 
cosmological quantities p x and w x (and hence the normalized Hubble expansion rate 
H) have reasonably simple forms which we can fit with the observational data. 

Let us first take into account the cases s = — 3n or s = for either of which Eq. (29) 
is linear in y 2 . However, for s = — 3n, Eq. (29) gives y oc y x , i.e., y = y + y 1 oc a~ 3n , 
which can never change sign for any finite value of a. Therefore we rule out s = —3n as 
a possible choice for a PDL crossing DE model. For s = 0, on the other hand, we obtain 
from Eqs. (26) and (29) that y = y 1 /(Cy 1 — 1) = constant, i.e., the corresponding fields 
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4> 1 and 2 are non-dynamical, and PDL crossing never happens either. Hence we rule 
out the choice s = as well. 

We next consider the case s = 3n/2 for which Eq. (29) is quadratic in y . 
Corresponding to < n < 1, we have < s < 3. The solution of Eq. (29) is given by 



Vo = a () a s 



i±,/i + ^ 



(30) 



where we denote a = l/(2Ca 1 ). From Eqs. (26) and (30) we see once again that the 
expression for y = y + y 1 cannot change sign, i.e., the PDL crossing is not possible. 

We are therefore left with the choices of the parameter s which render Eq. (29) 
to be cubic, quartic or even higher order algebraic equation in y %. The solutions for 
y are however very complicated and leads to messy expressions for the corresponding 
energy density. An inspection of Eq. (29) also reveals that the choices of s for which 
Eq. (29) is even order (quadratic, quartic, etc.) in y , can be ruled out anyway, because 
y — y + y 1 = Cy 1 y^ lrS / Zn can in no way change sign when the power (1 + s/3n) is an 
even integer. In a separate work which is in progress [50] we examine the cases for which 
Eq. (29) is cubic or higher order in y () . In the present paper, however, we choose to 
work with a simpler alternative by introducing one more k-field, in a way that Eq. (29) 
remains quadratic in y , and the above solution (30) could be used. 

4-2. Three k-field configuration 

Let us now consider the following configuration that involves three scalar fields: 

7 = 0,1,2 ; $ = (0 O ,0 1 ,0 2 ) ; rj (a) = -rjA a ) = v(a) , V 2 = • (31) 

From Eq. (23) we have 

*(l + JL) = yl(l-JL) = yl = _ 3". (32 ) 

The solution for y 2 is 

y, 2 = a 2 a~ 3n , a 2 = constant . (33) 

whence from Eqs. (19) and (20) we find the energy density and the pressure 
corresponding to the field (f> 2 as 

_ a -3(n+D + A2 ? p 2 = _^L. a -3(n+l) _ ^ ? (34) 



n + 1 n+1 
where A 2 is a constant. 

Once again, we assume the power-law solution (26) for y 11 so that the energy density 
and the pressure corresponding to the field <p 1 are given by (27), and y is related to y 1 
through the relation (29). Similar to the case of two k- fields, it is easy to show that the 
PDL crossing never happens for three k-fields as well, if the parameter s = —3n or 

§ There are also the choices s = — 3n/2 and s = — 6n for which Eq. (29) is quadratic in ^/y^ and y 
respectively. However, the solution for y corresponding to these choices is not qualitatively different 
from that corresponding to s — 3n shown above. The PDL crossing can never take place either. 
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so that Eq. (29) is linear in y . We therefore resort to the choice s = 3n for which Eq. 
(29) is quadratic in y and the solution for y is given by Eq. (30). Using Eqs. (26), 
(33) and (30) we obtain 



y = y +yi+V2 



i + 



2a 1 a 2s 



(35) 



(a + a 2 ) a s + o^a 5 ± a a~ 

a c ) = 0, which implies that the parameters a ,a 1 



For a PDL crossing at a = a c , y(a 
and a 2 should satisfy the relation 

(a 2 +a 1 a 2s ) 2 + 2a a 2 = 0. (36) 

Thus a smooth crossing from a w x > — 1 regime to a w x < — 1 regime is indeed possible 
for the choice s = 3n, provided one of the parameters a and a 2 must be negative with 
respect to the other, i.e., a a 2 < 0. However, a negative value of a 2 implies an eternally 
negative definite kinetic energy density p 2 for the field <f> 2 . A negative value a , on the 
other hand, does not necessarily imply a negative energy density for the field <j) , which 
is derived below. We therefore choose to work with a < 0, a 2 > 0. 

Let us now set s = 1/2, so that the DE perturbations are slowly propagating with 



speed c y 



n 



1/V6, and for which Eq. (20) could be integrated easily. Using 



Eqs. (19) and (20) one then obtains the expressions for the energy density and pressure 
components due to the field </> as 
6<x 



Po = 



7a 7 / 2 
7a 7 / 2 



1 ± 



, A(3a 8(3 2 a 2 
1 + + 



1±U + 



5 

6/3a 



15 



, 4,3a 
1 + + 



1 - (3af 2 
8p 2 a 2 



(l-[3a) 1/2 



(37) 



3 3 

where A is a constant and we denote (3 = —2a 1 /a . Using Eqs. (27), (34) and (37) 
we write the expressions for the total dark energy density (p x = p + p 1 + p 2 ) and total 
pressure (p x = P + P 1 + P 2 ), for the choice s = 3n = 1/2, as 

6a, 6a n 



Q(a +a 2/ 



7a 7 / 2 



+ 



± 



Px = 



5a 5 / 2 7a 7 / 2 
a 1 , a 



ABa 88 2 a 2 

1 + — + — 

5 15 



± 



1 + 



6/3a 



A(3a 
1 + + 



1 - fiaf 2 + A , 
8(3 2 a 2 



l-(3a) 1/2 -\, (38) 



7a 7 / 2 5a 5 /2 ^ 7 a 7/2 

where A = A + A 1 + A 2 . 

Now, from Eq. (6) we express the normalized Hubble rate (if = H/Hq) as 



H 



n 0m | 6(A + A 2 ) | 6A, 



± 



6A, 
7a 7 / 2 



1 



Af3a 



l-(3a) 3/2 +A, (39) 



a-' TaU 2 ' haF 2 

where A = a /p 0c ,A 1 = aJp 0c ,A 2 = a 2 /p 0c ,(3 = -2aja = -2AJA and A = \/p 0c . 
Eliminating A using the condition that at the present epoch t — t , a — 1 and H — 1, 
we obtain 



H 1 



± 



1 + O 0m (a 
6A 



1 



6(A +A a ) ^_ 7/2 



7/2 



1 + + 



7 
15 



(«- 7/2 " 1) 
1 - (3af 2 



6A 



-5/2 



1) 



, Ap 8/3 2 
1 + — + — 
5 15 



3/2 



• (40) 
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One may observe that the parameter (3 should not be greater than 1 in order that 
H 2 is real, and the model to be valid. Moreover, for positive fractional values of /3, 
for which the PDL crossing could take place, the range of validity of the model is up 
to a finite time in the future at which the scale factor is a max = 1/(3. Therefore, the 
lesser the positive fractional value of /3, the more the model could be extended to future 
epochs. 

In what follows, we fit this model with the observational data so as to determine 
the parameters A , A 1 and A 2 , and finally reconstruct the Lagrangian components 
P (y ), P-tiVi) an d P2i.ll 2) using the best fit values of the parameters. 



5. Fitting the model with observational data 

We use the most recent SNIa data released by the SN Cosmology Project (SCP) team, 
viz., the Union data-set [7]. It consists of a total of 414 SN samples analyzed from 
various sources, reduced to 307 most reliable data points after different outlier rejection 
cuts and selection tests. These 307 Union data includes large samples of SNIa from 
the older data-sets by Perlmutter et al, Tonry et al, Riess et al, and others [1, 2, 3, 4], 
SN Legacy Survey (SNLS) [5], as well as the distant SN observed by the Hubble Space 
Telescope (HST). The range of data extends up to red-shift z ~ 1.7, and the full data 
can be obtained at http:/ / supernova.lbl.gov/Union. 

These data provide the extinction-corrected distance modulus, given by 

Hobs(zi) = m obs (zi) - M , (41) 

where m b s is the apparent magnitude of the source SN and M is the absolute magnitude. 
Theoretically, the apparent magnitude is expressed as 

m th ( Zi ) = 5\og 10 [D L (z i )]+ M(M,H ), (42) 

where D^[zi) is the Hubble free luminosity distance, which for a spatially flat universe is 
defined in terms of the present matter density parameter Q 0m and the set of theoretical 
model parameters A = (A , . . . , A N _ X ) as 

D L (z t ) = (1 + Zi ) / ' (43) 

Jo tl [Zf, iZ 0m , AJ 

M is the magnitude zero point offset given by 



M(M, H ) = M + 5 log 10 



+ 25 = M-51og 10 /i + 42.38, (44) 



Mpc 

where h is the Hubble constant H in units of 100 Km s~ x Mpc -1 . 
The theoretical distance modulus fith( z i) is defined by 

Hth(zi) = m th (zi)-M = 51og 10 [D L {zi)\ +fi , (45) 

where 

fio = 42.38 - 51og 10 /i, (46) 
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is a nuisance parameter, independent of the data points, and has to be uniformly 
marginalized over (i.e., integrated out). 

The likelihood of flom and the theoretical model parameters A = (A , . . . , A N _J 
are determined by minimizing the quantity 

i=l i {Zi) 

where (Ti(zi), (i = 1, . . . ,N tat ), are the total uncertainties on the data points. For the 
Union data-set [7] N tot = 307. 

Now, the marginalization with respect to /io can be done by the following two-step 
procedure [51]: 

1. Expanding x 2 in terms of /io as 

X 2 (/i ;fi 0m ,A) = A(n 0rn ,A) - 2/2 B(n 0m ,A) + ulC, (48) 

where 

a, a \ \j*obs(zj) - Ufa no = 0)] 2 

i=l a i ^ Zi > 

Ntot i \ / _ n\ Ntot i 

B(f!o - A) = h sfw ' c= lr?F(i)' (49) 

one finds that Eq. (48) is minimum for /i = B/C. 

2. Substituting this value of m Eq. (48) we obtain 

X 2 (n 0m ,A) = x 2 (^o = ^om^ = A(n 0m ,A) - i]3 2 (Q 0m ,A). (50) 

Now instead of minimizing x 2 in Eq. (47), one can minimize x 2 in Eq. (50) with respect 
to Q 0m and A = (A , . . . , A N _A, since both the minimized values would be the same. 

In the case of the three k-field (N = 3) model discussed above, corresponding to 
the choice s = 3n = 1/2, we find good fits to the SNIa data for positive fractional values 
of /3 = —2A 1 /A () , and with the negative sign of square roots appearing in the expression 
for H 2 in Eq. (40). We fit the observational data for two parametric choices of f3 — 0.5 
and 0.25, for which the model remains valid up to sufficiently large future red-shifts 
z = —0.5 and —0.75 respectively. 

j3 = 0.5: With this choice we get the minimized value of x 2 = 310.411 for the 307 Union 
data-set [7], and the best fit values Q 0m = 0.26662, A 1 = 1.46422, A 2 = 0.159546. The 
minimized x 2 is better by A% 2 = —2.69 from that found with the Union data for the 
two-parameter flat ACDM model [7, 8]. The minimized x 2 P er degree of freedom (dof) 
is also improved by A(x 2 /dof) = —0.0019 from that found in ref. [6] by fitting the flat 
ACDM model with the older SNIa gold data-set [4]. 



j3 = 0.25: With this choice we get the minimized x 2 = 310.385, which is improved 
further by Ax 2 = —0.26 than that obtained for the previous choice {(3 = 0.5). The 
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best fit values of the present matter density parameter and the model parameters are 
obtained as: tt 0m = 0.26103, A 1 = 3.58355, A 2 = 0.174564. 




Figure 1. Evolution of w x (z) and p x (z) (best fit with the Union data, plus 
la errors) for < z < 1.75, as well as the la (white), 2a (light grey) and 
3a (dark grey) contours in A l — A 2 parameter space, for the choices j3 = 0.5 
(upper panels) and (3 = 0.25 (lower panels). The point z c denotes the value 
of the red-shift at which the PDL (w x = —1) is crossed. The best fit points 
for both the choices are shown by the dots at the middle of the la A 1 — A 2 
contours. The cosmological constant, which corresponds to A 1 = A 2 = 0, is 
found to lie on edge of la A 1 — A 2 contour for both the choices. 

The evolution of w x (z) and p x (z) (in units of p 0c ), alongwith the corresponding la 
errors, for < z < 1.75, as well as the la, la and 3a A 1 — A 2 contours are shown in 
Fig. 1, for both the choices (3 = 0.5 and 0.25. One may also have the Qo m — A l and 
Qo m — A 2 contour plots for either choice, however for brevity we have suppressed those 
plots. The cases (3 = 0.5 and (3 = 0.25 are potentially not much distinct. The w x (z) 
vs z plots show that the value z = z c at which the best fit w x crosses —1 is 0.2125 for 
f3 = 0.5 and 0.1988 for (3 = 0.25. Both these values of z c are fairly in agreement with 
those obtained for popular model-independent ansatze for w x (z) [13], H(z) [14], p x (z) 
[15], q(z) [16], etc. The cosmological constant (corresponding to A 1 = A 2 = 0), is found 
to be about la away from the best fit result for both the choices (3 = 0.5 and f3 = 0.25. 

The left panel of Fig. 2 shows the variations of the DE EoS parameter w x (z), the 
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Figure 2. The left panel shows the extrapolations of w x (z),w(z) and q (best 
fit with the Union data) to the red-shift range —0.5 < z < 2.5. The DE 
EoS parameter w x shows a tendency to back from a minimum value in future 
towards the cosmological constant value —1. The right panel shows a nearly 
tracking nature of the dark energy density p x (extrapolated best fit) relative 
to the matter density p m for high red-shifts z > 1. Both the plots correspond 
to the choice (5 = 0.25. 



total EoS parameter w(z), Eq. (10), and the deceleration parameter q(z), Eq. (11), for 
the best fit parameters flo m , A 1 ,A 2 , extrapolated to the red-shift range —0.5 < z < 2.5. 
The transition from the decelerating (q > 0) phase to the accelerating (q < 0) phase 
occurs at z ~ 0.45 and the values of both w and q are greater than —1 at present 
(z = 0). Whereas the present value of w x < —1, there is a tendency that w x would 
turn back from a minimum value towards the cosmological constant value (= —1) in 
distant future (z < —0.5). Thus the phantom (w x < —1) phase may be transient, 
so that future singularities, such as the Big Rip, may be avoided. The right panel of 
Fig. 2 shows the variations of the matter density p m and the extrapolated best fit DE 
density p x (both in units of the present critical density poc) with z, for —0.7 < z < 2.5. 
The DE density is found to almost track the matter density for higher red-shifts, until 
equalling it at z ~ 0.4 and dominant thereafter ||. Thus the model apparently resolves 
the coincidence problem, side by side providing a PDL crossing DE. The extrapolation 
of the best fit p x to future epochs also show that the DE density, though dominant, 
does not run away to very high values, causing Big Rip-like singularities. Both the plots 
of Fig. 2 correspond to the choice (3 = 0.25, for which the model fits better with the 

|| For very high red-shifts (z > 10) we find that extrapolated p x exceeds p m once again and the DE 
becomes more and more dominant as we go further and further in the past. However, the best fit model 
parameters, which give rise to such a variation of p x (z), are found by fitting the Union data available 
only up to z ~ 1.7. Therefore the extrapolation of p x (z) to such high red-shifts, using these parametric 
values, cannot be trusted with much confidence after all. 
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Figure 3. The upper left panel shows the evolution of the individual EoS 
parameters w (z),w 1 (z),w 2 (z) (extrapolated best fit with the Union data) 
for —0.5 < z < 2. Whereas w 1 and w 2 are negative and positive constants 
respectively, w changes from a negative to a positive value in very recent 
past. The upper middle panel shows the evolution of the individual energy 
densities p (z), p 1 (z), p 2 (z) (extrapolated best fit with the Union data) for 
— 0.7 < z < 1.5. While p 1 and p 2 are positive definite, p changes from a 
negative to a positive value at a distant future epoch (z « —0.5). The upper 
right panel is a parametric plot (for < z < 1.75) showing the variation 
of the k-essence Lagrangian p x = P (reconstructed using the best fit model 
parameters) with y = y + y 1 +y 2 - The lower panels show the parametric plots 
(for < z < 1.75) of the reconstructed normalized Lagrangian components 
P ni = Pj/PW versus the normalized field gradients y ni = y T /y^ for I = 0, 1, 2. 
All the plots correspond to the choice (3 = 0.25. 



Union data, as mentioned above. 

All the plots of Fig. 3 also correspond to the choice (3 = 0.25. The upper left panel 
of Fig. 3 shows the extrapolated best fit variations of the EoS parameters Wj = PJ p t 
corresponding to the individual fields f (/ = 0, 1, 2) with red-shift z, for —0.5 < z < 2.5. 
While w 1 = —n = — 1/6 and w 2 = n = 1/6 are negative and positive constants equal 
in magnitude, the EoS parameter w increases from a negative constant value (equal to 
w 1 = —1/6) in distant past, crosses zero in very recent past (z ~ 0.05) and then goes on 
increasing with positive values. The upper middle panel of Fig. 3 shows the extrapolated 
best fit variations of the individual field energy densities p 7 with z, for —0.7 < z < 1.5. 
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p 1 and p. 2 are positive definite and decreases with time (i.e., with decreasing z). p is 
negative in the past and at present, but increases with time, and may become positive 
in distant future (z < —0.5), as the extrapolations show. In magnitude the value of p 2 
is found to be very small compared to p and p 1 for all the red-shift range shown in the 
figure. p and p i; on the other hand, are nearly equal in magnitude for < z < 0.5. 
Now, the effective DE EoS parameter w x is given by 

w p + w 1 p 1 + w 2 p 2 , . 

^ = — — — , P X =P +P 1 +P2- (51) 

Px 

Since p x is always positive, the negative definite and positive definite contributors to 
w x are (p 1 (for which w 1 < 0,p 1 > 0) and <p 2 (for which w 2 > 0,p 2 > 0) respectively. 
The field O , on the other hand, gives positive contribution (w p > 0) to w x right from 
the past up to z ~ 0.05 and negative contribution {w p < 0) thereafter. However this 
field O , though has a negative energy density (phantom-like), could not have led to the 
PDL crossing at z ~ 0.2, as its negative contribution to w x begins only after z ~ 0.05. 
It is therefore the (non-phantom-like) field (f> 1 , with positive definite energy density p i; 
which is effectively responsible for the PDL crossing at z ~ 0.2. This is in contrast with 
the usual quintom mechanism [25] where the negative energy density (phantom) field 
leads to the PDL crossing at late times. Moreover, none of the three fields above has 
EoS parameter w T less than —1, as in quintom models. 

We reconstruct the multiple kinetic k-essence Lagrangian P = p x using the best 
fit values of the parameters f2 0m , A 1 , A 2 for the choice f3 = 0.25. The parametric plot of 
p x against the total temporal field gradient y = y +y 1 +y 2 , for red-shifts z = — 1.75, 
in the upper right panel of Fig. 3 shows three separate phases: (i) p x > for y > 0.18, 
whence the DE is decelerating (w x > 0), (ii) —0.7 < p x < for < y < 0.18, whence 
the DE is accelerating but with w x > —1, (iii) p x < —0.7 for y < 0, whence we have 
a super-accelerating (i.e., phantom-like) DE with w x < —1. The parametric plots of 
the reconstructed normalized k-essence Lagrangian components P ni = PJP^ (where 
p(o) _ | z=0) / — 0, 1, 2) against the corresponding normalized temporal field gradients 
Vn! = yjyf^ (where yf^ = yj\ z =o, I = 0, 1, 2), for < z < 1.75, are shown in the lower 
panels of Fig. 2. In the reconstructions, we have set for simplicity A x and A2 appearing 
in Eqs. (27) and (34) to be zero, whence the cosmological constant term in Eq. (38) 
is given by A = A . The pressure P ni is found to fall off with increasing values of 
y m , whereas the pressures P no and P n2 are found to increase, but with opposite rates, 
with y no and y n2 respectively. Since the best fit values of both the parameters A 1 , A 2 
are positive, the only parameter with a negative value is A = —2A 1 /j3 = —8A 1 (for 
(3 = 0.25). Therefore, the fields which contribute to the negative part of w x are <pi (due 
to the negative constant value of EoS parameter w 1 , but positive definite energy density 
pj and 0o (due to the negative parameter A , or a , which drives y = O and/or P 
to be negative, in some regime). However, as explained above, the field O alone is not 
responsible for the PDL crossing, and it is the field <p l with Lagrangian component P 1 
which effectively causes the super-acceleration {w x < — 1). 
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6. Conclusions 

While a PDL crossing DE model involving a single k-essence scalar field suffers from 
cosmological instabilities [35] , a set of multiple k-essence fields are shown to give rise to 
such PDL crossing, satisfying the conditions for stability even in the simplistic scenario 
where the k-fields are mutually non-interacting and the individual k-essence Lagrangian 
densities are purely functions of the kinetic terms. Although there are severe restrictions 
imposed on the k- field configuration by virtue of the stability criterion, i.e., c 2 x > 0, and 
the non-super-luminescence (c x < 1) of the DE density perturbations, we still find 
certain field configurations for which the DE EoS parameter can cross the cosmological 
constant barrier. We have worked out the general equation of constraint on the field 
configuration, with an appropriate choice of the propagation speed, viz., c 2 = a constant, 
lying between and 1. While a configuration of two kinetic k-fields, satisfying the 
constraint equation, leads to complicated solutions, we have shown that a three k-field 
configuration gives rise to a cosmologically stable PDL crossing DE model. For certain 
choices of some of the parameters, the model is found to fit well with the latest SNIa 
Union data compiled by the Supernova Cosmology Project team [7]. The value of the 
red-shift z = z c at which the DE equation of state parameter w x (best fit with the 
data) crosses the PDL is found to be ~ 0.2, whereas the best fit present matter density 
parameter is found to be Qo m ~ 0.265. These are in good agreement with the values of 
z c and Q 0m found in the likelihood analysis done using other model-dependent or model- 
independent ansatze in the literature [16, 17, 49]. On the other hand, the value of z at 
which the universe transits from a decelerating phase to an accelerating phase is found 
to be ~ 0.65, which also agrees fairly well with other independent studies [16, 17, 49]. 

In the model which we have presented here, the value of the square of the speed 
of propagation of the DE density perturbations, c 2 x , is as low as 1/6 throughout 
the entire course of the DE evolution. We have worked out the expression for c 2 x 
in a general setting of multiple k-essence fields considering the most general scalar 
type of metric perturbations in a Friedmann-Robertson- Walker background. Although 
the perturbation equations have been extremely coupled in the curvature and the 
isocurvature modes, assuming the coupling to be very weak in the sub-Hubble scales the 
speed of propagation of the curvature perturbations in the multi-k-field case is found to 
be the one that generalizes the well-known expression for c x for a single k-essence field 
found in ref. [45]. 

The present model exhibits some additional features: Firstly, the dark energy 
density (best fit with the Union data) almost tracks the matter density in distant 
past until exceeding the later at z ~ 0.4, thus providing a possible resolution to the 
coincidence problem. Secondly, the best fit dark energy density, extrapolated to future 
epochs, does not run away to very high values leading to Big Rip-like singularities even 
in distant future. Thirdly, the extrapolated best fit of the dark energy equation of 
state parameter shows a tendency to turn back from a minimum negative value (< — 1) 
towards the cosmological constant value (= —1) in distant future, thus indicating a 
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transient phantom regime. 

We have reconstructed the entire k-essence Lagrangian P(y ,y 1 ,y 2 ), as well as the 
individual Lagrangian components P I {y I ') in terms of the corresponding field gradients 
y l = <pj, using the values of the model parameters best fitting the observational data. 
The reconstructed total Lagrangian, P is found to be a smooth function of the total field 
gradient y — YlVn i- e -> t nere are no discontinuities or multi-valuedness in the parametric 
plot of P versus y for the red-shift range under consideration, < z < 1.75. All the 
individual Lagrangian components P t as functions of the corresponding y t also show 
the same smooth nature. The reconstructed Lagrangian P is therefore substantially 
better than the double-valued Lagrangian found for a single kinetic k-essence field in 
ref. [52] using popular parameterizations, such as that due to Alam et al [14] or due to 
Chevallier-Polarski-Linder [13] and with the old SNIa gold data-set [3]. 

A few problems that may be looked upon in the context of the present model are 
as follows: (i) Can we find a systematic prescription in order to build a stable PDL 
crossing multiple kinetic k-essence model, instead of choosing the number of k-fields 
by hand and implementing the value of c 2 x arbitrarily? (ii) Can we extend the model 
to sufficiently higher red-shifts so as to match with the CMB and SDSS data? (iii) Is 
there a way to ascertain how far the isocurvature modes of the density perturbations 
affect the model at higher red-shifts? (iv) Can we ascertain the status of the future Big 
Rip singularities, if any, in the context of multi-k-essence, generically, i.e., by not just 
resorting to particular models? Or, can we generically ascertain whether the phantom 
regime would be eternal or transient for multiple k-essence? Research works addressing 
some of these questions are already in progress [46] and we hope to report on them in 
near future. 

Finally, it should be mentioned that although we have not specifically looked 
to construct tracking/scaling dark energy models in the present paper, the possible 
avoidance of the coincidence problem in our multi-k-essence model is particularly 
noteworthy. The value of c x remains constant and less than the speed of light, unlike in 
the coincidence resolving single field k-essence models, which are shown to be affected 
by the super-luminal propagation (c x > 1) of the density perturbations [47]. However, 
whether our purely kinetic model specifically belongs to the class of tracking or scaling 
k-essence dark energy, is an interesting issue which we hope to circumvent in future. 
Assisted accelerated solutions are shown to exist generically for the multi-field k-essence 
models that have scaling solutions [53]. The nature of such assisted acceleration may 
also be of some interest for the purely kinetic multi-k-essence, that we have studied here, 
and is worth investigating. 
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Appendix: Speed of propagation of multiple k-essence cosmological 
perturbations 

Let us consider the system of N mutually non-interacting k-essence fields coupled to 
gravity. The gravitational field equations for such a system are 

N-l 

G^u = k 2 T^ u , T^y = T^ v , (52) 

1=0 

where TL,, (I = 0, . . . , N — 1) are the energy-momentum tensor components due to the 
individual k-fields <& = (0 O , • • • , Ar _ 1 ), with kinetic terms X = (X , • • • , X N _ 1 ), where 

The line element with the most general scalar type of perturbations, in an arbitrary 
gauge, is given by [54]: 

ds 2 = - [l + 2A(x\t)] dt 2 + 2a(t) diB(x\t) dx l dt 

+ a 2 (t) [(l + 2i/j(x\ t)) 8ij + VjdiEix*, t)] dx*dx j , (53) 

where a(t) is the cosmological scale factor, and A(x l , t), B(x\ t),tp(x\ t) and E(x\ t) are 
the perturbed order variables. The perturbed energy-momentum tensor components are 
given by 

T° = -(p + 5p), T^~(p + p)d iV: 

Tl = (p + 5p) 6> + 1 (V, V" + ^ vr , (54) 

where p and p are the total background energy density and isotropic pressure 
respectively, p = p + 5p and p = p + 5p are the perturbed total energy density 
and isotropic pressure respectively, v is the velocity (or, flux) related variable, n is 
the anisotropic stress, and k is the co-moving wave-number. Decomposing into the 
components due to individual k-fields, we have 

iii i 

p = J2 p i - p = J2 p i - (p+p) v = J2( p i +p i) v i ' 7r = $^ 7r i- ( 55 ) 
ii i i 

The background field equations are the Friedmann equations and the equations of 
continuity of the individual scalar fields: 

H 2 - f-V = ^P, (56) 



a 2 k 2 



H = - 2 = - -(P + P) , (57) 

a a A 2 

Pj + 3i7(p J + P / ) =0 , for all I = 0,...,N -1, (58) 
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The equations for the cosmological perturbations 

SGi = k 2 5T^ (59) 

can be decomposed as 
k 2 k 2 

if £ — —ip = — — 5p : 5Gn component , (60) 
a 2 2 

£ - -^x = ~Y %(P + P) V :SG i com P onent ? ( 61 ) 

X + H X = A + ip :5G{- ^5j5G k k component , (62) 

£ + 2H£ + (3H - — ) A = — {5p + 35p) :5G k k -5G° component, (63) 



A + V • (66) 



where 

7 2 

X = a(B + aE) , £ = 3 (HA- j>) + ^ X, (64) 

and note that tt = for a scalar field system in RW background. The quantities ip, x 
and £ correspond respectively to the three-space curvature, the shear, and the perturbed 
expansion of the normal-frame vector field [39, 40]. Eq. (60) gives the Arnowitt-Deser- 
Misner (ADM) energy constraint, Eq. (61) gives the momentum constraint, Eq. (62) 
is the trace-free part of the ADM propagation, and Eq. (63) is the Raychaudhuri 
equation [40]. In addition to the above equations, we have the energy and momentum 
conservation equations given by 

5 Pl + 3H (5 Pl +6P I )= ^ - 3HA - K)j ( Pl + Pj , (65) 

frlVfo + ^K] k f 6P X 

^(fij + Pj) a\ Pl + P^ 

The entropic perturbations are given by 

e, = SPj - c 2 i 5 Pl , e = 5p - c 2 s 5 P , (67) 

where c 2 = PJ Pl and c 2 = p/p. 

Under a coordinate transformation 

x a -> x a + q a , q a = (q°, q l ) ; q° = St , q l = 5 ij dj5x , (68) 
the metric perturbations transform as 

A^ A-5i , B ^ B- — + a5x , ip -> i/j - H6t , E -> E - 5x , (69) 

a 

whereas the field perturbations transform as 

5(f); -> 50, - j> z 8t . (70) 

Since we are always free to choose St and 5x, we can impose two gauge conditions on the 
functions A, B, ip, E or Sp or diV. We choose to work in the commonly used longitudinal 
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gauge in which B = E = 0. Then from Eq. (64) x — 0, whence Eq. (62) gives A = —if), 
and the expression for £ in Eq. (64) reduces to 

£ = - + ■ (71) 

Curvature Perturbations : 

The so-called gauge-invariant co-moving curvature perturbation variable is given 
by [55] 



<E> = ip — -Hv . 

k 



(72) 



From Eqs. (60), (61), (67) and (71), one derives the equations for the curvature 
perturbations as 



ip + Hip 
H 2 



(p + p)at \H 



k 2 a 
Y k 



5p v 



(p + p)v , 



k 2 



r<*>, 



where the gauge invariant variable 5p v is defined by 



d~Pv = 6p- -pv = - — 
Isocurvature Perturbations: 



1 ( P + P k 



H 



$ + e . 



(73) 
(74) 

(75) 



The total entropic perturbation can be decomposed as [40, 41] 



, 2W 1 (Pi + Pi) (Pj + Pj) „2 o 

/ i 1 i,j p + p 

„2 _ „2 „2 



(76) 



where c 2 7 = c 2 — c 2 , and the so-called gauge-invariant isocurvature perturbation variable 

Spj 5pj 



Sjj is defined as 



" Pj + Pi Pj + Pj 



(77) 



From Eqs. (65), (66) and (67) one obtains the equations for the isocurvature 
perturbations 



S rJ = 



k 



a 



V IJ - 3 # e /.7 , 



V IJ = V J- V J , e IJ 



Pj + Pj Pj + Pj 



(78) 



S„ + 21fS„ - ■ ™ 



Aj 



Sjj + E^^(^+^) 



+ 



2a 2 



<jSjj + cLE 



K 



K 

Pk + P K 
P + P 



- 3He ZJ - 3 (H + 2# 2 ) e„ + 

2 hv 



P + P 
(S IK + S JK ) 
9H 2 



c 2 e 



A' 



P + P 



( e iK + ^Jk) 



k 2 

a 



2 \ °IJ 



P + P 



+ e 



u i ' 



(79) 
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where c 2 = c 2 + c 2 . 

u i j 

Eqs. (73), (74) and (78), (79) form an extremely coupled set of differential equations 
in the curvature and isocurvature modes, very difficult to solve in general for any system 
of multiple fields. However, leaving aside the isocurvature terms, one may still work out 
the effective speed of propagation of curvature perturbations in such a system. In fact, 
in the large scale limit (k/a — > 0), the curvature mode couples very weakly to the 
isocurvature ones, as is evident from the term (k 2 /a 2 )c IJ 8p v /(p +p) on the right hand 
side of Eq. (79). For a detailed discussion on the issue of decoupling of isocurvature 
and curvature modes in multiple fluid or field systems, see ref. [40]. 

Now, in the case of mutually non-interacting multiple k-essence scalar fields, we 
have the expressions for the total energy density and pressure 

N-l N-l 



dP 

X^-P^X,) , (80) 

7=0 1=0 ' 

N-l 

1=0 



Using the continuity equations for the fields $ = (0 O , • • • , <fi N and the fact that 
Xj = 2 /2 in RW background, we obtain the expressions for the perturbed variables 5p 
and 5p as 



N-l 



5 P =J2 s Pi > s Pi = ^rh ( 3H % a * ' & < ~ ^ ' (82) 



/=0 
N-l 



3HS 2 



5 P =J2 8P, , 5P i = {fj ( 3Hc >i - & i - ^) » ( 83 ) 

7=0 

where we have defined a 1 = and c 2 = ^-/§x~- 

The expressions for v r and v, obtained from the — i components of the energy- 
momentum tensors, are given by 

k k ^ 

= v = - y,p>°>- ( 84 ) 

' 1=0 

Defining further a = °~i an d °~u — °~i — °~. n we get from Eqs. (72), (73) and (75) 
$ = V " ^ + 0(a XJ ) , (85) 

i> + = ^ + 0(a u ) , (86) 

6n 1=0 °i 

where N is the total number of fields and for the reasons mentioned earlier we have left 
aside the terms depending only on a u (and/or its time derivatives), which give rise to 
the isocurvature modes. 
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The equation (73) for the curvature mode can now be expressed as 



H 2 



(Px +Px) a3 



(_Px +Px) q3 



+ 



where 



P 



N-l . 



/=0 I 



N-l 

E 



p I +p I 



c 2 



(89) 



(Px + P> 

is the effective speed of propagation of the curvature mode of scalar perturbations 
[46]. For a single k- field (N = 1,1 = 0), Eq. (89) reduces to the well-known result 



X 

~.2 . 



C 2 

i 

■'I 



2 = 8P /dX [45], whereas for purely kinetic multiple k-essence fields, we have 
<r = Pji ' p I , whence c 2 = pj XXPj/A)- ^ an Y °f ^ ne fields is non-propagating, i.e., 



z 2 = for some value of /, then c 2 x 
Under the transformation 



u = C$ 



c = 



a VPx + P> 

c v H 



the above equation (88) reduces to the familiar form [45, 39, 54, 55] 



u" + 



2 i 2 C 



(90) 



(91) 



where the prime {'} here denotes differentiation with respect to the conformal time 77 
(=fdt/a(t)). 
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